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Abstract 

We analyze and simulate a two dimensional Brownian multi-type particle sys- 
tem with death and branching (birth) depending on the position of particles of 
different types. The system is confined in the two dimensional box, whose bound- 
aries act as the sink of Brownian particles. The branching rate matches the death 
rate so that the total number of particles is kept constant. In the case of m types 
of particles in the rectangular box of size a, b and elongated shape a > 6 we ob- 
serve that the stationary distribution of particles corresponds to the m-th Laplacian 
eigenfunction. For smaller elongations a > b we find a configurational transition to 
a new limiting distribution. The ratio a/b for which the transition occurs is related 
to the value of the m-th eigenvalue of the Laplacian with rectangular boundaries. 
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I Introduction 

It is remarkable how the simple systems with few deterministic rules (such as 
Life [1] or cellular automata [2]) can generate complex structures. In the population 
dynamics however one often uses stochastic models. For example as has been re- 
cently shown [3] the addition of stochastic factor into the Life game favors diversity 
of structures, contrary to the original model in which diversity is the decreasing 
function of time. Introduction of the probabilistic factor in the cellular automata 
in the description of the dynamics of the social impact in the population [4] leads 
to the complex spatial and temporal intermittent behavior. In the genome popu- 
lation dynamics [5-6] one uses stochastic processes such as super-Brownian motion 
or Fleming-Viot processes. The model presented in this paper is a special type of 
population dynamic stochastic model. 

The dynamics of systems with two competing species has been studied with 
the emphasis on the spatial heterogeneity influence on the temporal evolution and 
spatial organization [7-9]. In the case of strong competition only one of the two 
survives, which means that the average lifetime of two species can be different. In 
our model we are concerned with the spatial distribution of three or more competing 
species. Our model differs from those mentioned above because we study the case 
when the total number of particles is constant and the average lifetime of different 
species is the same in the long time limit. 

Here we describe a behavior of a population of m different types of particles 
with Brownian dynamics confined in the two dimensional box, whose walls act as 
sinks for the particles. Additionally we assume that if two particles of different type 
occupy the same lattice point, both are killed. The birth rules are chosen in such a 
way as to keep the number of particles constant at each time step and to ensure that 
the average lifetime of any type of particles is the same in the long time limit. As an 
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example we show (Fig.l) the stationary configurations for 3 types of particles in a 
rectangular box of size a x b. For a/b > 1.63 the particles of different types occupy 
domains of rectangular shapes (Fig. la). We call such configuration elementary 
and show (in section V) that it is related to the third Laplacian eigenfunction. 
When the size ratio decreases below 1.63 the configuration changes its character as 
shown in Figs.lbc. Here the domains have the shapes that are not related to the 
Laplacian eigenfunctions. We shall call this transition the configurational transition. 
In Section VII we show that the ratio of rectangle sides a/b at the transition can 
be obtained from the simple condition involving the third Laplacian eigenvalue. In 
a natural way our model provides the probabilistic interpretation of the Laplacian 
eigenfunctions. We would like to emphasize that our stochastic model leads to a 
deterministic limiting distribution, whereas for example in super-Brownian motion 
and Fleming- Viot processes the limiting distribution has the fractal nature. 

The paper is organized as follows. In Section II we briefly discuss the discrete 
analogues of super-Brownian motion and Fleming-Viot processes. In Sections III 
and IV we describe in detail our model. The connection between the stationary 
state of the model and the Laplacian eigenfunctions is given in Section V and com- 
puter simulations are described in Section VI. The analysis of the configurational 
transition and the concluding remarks are contained in Section VII. 

II Super-Brownian and Fleming-Viot processes: 
particle systems with death rate independent of position 

Super-Brownian motion and Fleming-Viot processes are usually discussed in 
the continuous time and space state setting. We will present their discrete analogues 
for the purpose of comparison with our own model introduced in Sections III and 
IV below. In the first model we consider particles on the two dimensional square 
lattice. At every time step t = 1, 2, 3 • • •, each particle either dies or branches off 
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into two offsprings with probability 1/2. If the particle branches, both offsprings 
occupy the same lattice site as the parent particle and then each goes to one of 
the four neighbour lattice sites with probability 1/4. The events are independent 
for all the particles in the population. Suppose that at time t = 1 the particle 
system consists of j particles and every particle is located at (0,0). Let X J S be a 
measure-valued process whose value at time s is defined as follows. The measure 
X 3 a (A) of an open subset A of R 2 is equal to the number of particles at time t = [s] 
which lie in y/JA ([s] is the integer part of s). Consider the sequence of processes 
{Xj u /j, u > 0}j>i where s = ju and u plays the role of the rescaled time. This 
sequence of processes converges as j — > oo to a measure-valued diffusion called 
super-Brownian motion with the initial state 5(o,o) (mass 1 concentrated at (0,0)). 
The limiting distribution of the process has the fractal nature in dimensions d > 2. 
At any fixed time, the state of super-Brownian motion is a measure whose support 
has Hausdorff dimension 2 [10], for d > 2. In other words, the volume occupied by 
the particles with the linear size of the system, L, scales as L 2 irrespective of d > 2 
(L — > 0). For d = 1, the limiting distribution of the process at a fixed time has a 
continuous density. 

The second model differs from the first one in that the population size is fixed 
and equal to j. The dynamics are now the following. First suppose that k = 
1,2, — 1 and n > 1. In order to obtain the state of the process at time 
t = nj + k + 1 from that at t = nj + k, we choose randomly one particle and kill 
it. Next, another particle is chosen from the surviving ones and it branches into 
two offspring which occupy the same lattice site as the parent particle. If t = nj 
then we obtain the new configuration at time t = nj + 1 by letting each of the 
particles move to one of the 4 nearest sites on the lattice, with probability 1/4, 
independent of all other particles. We renormalize the system in order to obtain 
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a continuous limit. Suppose that at time t = 1 the system consists of j particles 
located at (0, 0). Let X\ be a measure defined as in the first model, i.e., the measure 
Xj (A) of an open subset A of R 2 is equal to the number of particles at time t = [s] 
which lie in y/]A. Then the sequence of processes {X? 2u /j, u > 0}j>i (s = j 2 u) 
converges as j ' — > oo to a measure-valued diffusion called the Fleming-Viot process 
[6, 11] with the initial state 5(o,o)- This process has the same fractal nature as 
the super-Brownian motion. More precisely, the state of the Fleming-Viot process 
at a fixed time is a measure which has the same Hausdorff dimension (and other 
fractal properties like Hausdorff measure, packing measure, etc.) as super-Brownian 
motion. The Fleming-Viot process is the super-Brownian motion when the latter is 
conditioned to have a constant total number of particles. 

Recently there has been growing interest in models incorporating dependence of 
the motion of individual particles on the current configuration [11,12]. We propose 
to study a model with a constant population size in which particles can die and 
branch. In our model, the death of a particle will depend on its location and thus 
it differs from the two aforementioned models. The simplest case is when a particle 
dies if and only if it moves to a set of the designated sites on the lattice. 

Ill Particle system with death depending on position 

We fix a connected subset D s of the square lattice with the mesh size e, denoted 
(eZ) 2 . The particles in our model die if they move outside _D e , so D e plays the role 
of the state space. The number of particles is fixed and equal to j. Transitions from 
the state of the system at time t = k to that at time t = k + 1 may be described as 
follows. First each of the particles goes to one of the 4 nearest sites on the lattice 
(eZ) 2 , with probability 1/4, independent of all other particles. Then all particles 
which are outside D e die. An equal number of particles is chosen uniformly from 
among the surviving particles. Each of the chosen particles splits into two offspring 
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which occupy the same site as the parent particle. Hence, the number of particles 
in our model is constant between generations. 

Fix some open connected set D C R 2 and let D e = D n (eZ) 2 . Suppose that 
at time t = 1 each of the j particles occupy sites in D e . Let X J S ' £ be the measure 
valued process whose value at time s is defined as follows. The measure X J S ,S (A) of 
an open subset A of R 2 is equal to the number of particles which are in A at time 
[s). 

The qualitative long time behavior of our system is much different from that 
in the case of the super-Brownian motion or Fleming-Viot process. A typical par- 
ticle configuration in both models discussed in the previous section has a frac- 
tal nature. Rigorously speaking, the limiting continuous models are measure- 
valued diffusions whose states are measures supported on fractal sets [11]. In 
our model, increasing the number of particles j and decreasing the mesh e of the 
lattice results, in the long run, in a stable distribution which is a suitably nor- 
malized first eigenfunction of the Laplacian on D with zero boundary values. In 
other words, if f(x,y) denotes the first eigenfunction of the Laplacian with zero 
boundary values in D then hmX^ E2 (dx,dy)/j = cf(x,y)dxdy or, more precisely, 
\\mX 3 s ^ e2 (A) / j = f A cf(x, y)dxdy for every open set A C R 2 , where < c < oo and 
the limit is taken as e — > 0, j — > oo and s — > oo. Here c = 1/ f D f(x, y)dxdy. 

One physical interpretation of the first eigenfunction is that it represents the 
probability distribution, after a long time delay, for a randomly moving particle 
conditioned to stay within the domain [13]. The normalization of the eigenfunction 
is necessary to make the total probability equal to 1 (in the case of probabilistic 
interpretation) or to make its integral equal to the total mass of particles in our 
model (we normalize the mass by dividing the measure X 3 ^ by j). 

We offer a heuristic argument showing the convergence of distributions in our 
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model. The remarks are not meant to be a rigorous proof — that does not seem 
to be trivial and will be the subject of a forthcoming paper [14]. Notice first that 
because of the diffusive scaling x — > ex and s — * e~ 2 s, each particle, in the limit 
e — > 0, executes a Brownian motion in D with a jump, upon exiting £), to a point 
occupied by a fellow particle chosen uniformly at random. Second, since particles 
interact only through the boundary of D by a random choice from the remaining 
particles, the equal time pair correlations are inversely proportional to the total 
particle number j, and therefore the particles are uncorrelated in the limit j ' — > oo. 
Thus the limiting measure X s = lim e jX^ J e2 /j exists and is deterministic. Let us 
express this limit via its density X S (A) = J A £(s; x, y)dxdy. Since all particles reside 
in D it follows that £(s; x, y) > and it vanishes for points on the boundary of D. 
Now the average exit time of a typical particle from D is the reciprocal of Ai, the 
first eigenvalue of the Laplacian in D with zero boundary values [13], which shows 
that the per particle rate at which jumps take place is exactly Ai. Thus the density 
£(s; x, y) is the solution of a heat flow problem in D with a heat source of strength 
\i£(s;x,y) and absorbtion at the boundary, i.e., 

dt/ds + A^X^. (3.1) 

Here, A = —l/2(d 2 /dx 2 + d 2 /dy 2 ) is the Laplacian. As s — > oo, the density 
converges to a solution of the stationary problem for which the normalized first 
eigenfunction is the unique non- negative solution having total integral 1. 

IV Multi-type particle system 
The first eigenfunction of the Laplacian with zero boundary values already has 
a natural probabilistic interpretation [13] and the model described in the previous 
section provides a new one. It seems that so far the higher eigenfunctions do not 
have a natural probabilistic interpretation. A model described in this section may 
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be a first step towards such an interpretation. 

Fix a connected subset D e of the square lattice (sZ) 2 with the mesh size e. 
In this model, each particle will reside in D £ and have one of m possible types 
k = 1, 2, . . . , m. Typically, at each time t = I some particles will be chosen to split 
into two offspring. In such a case we will say that a new offspring was born at time 
t = I and if this particle is killed at some later time t = n then we will say that the 
lifetime of this particle was T = n — I. The transition mechanism of the system, 
which depends on the positions, types, and lifetimes of the particles, is the following 
one. First each of the particles goes to one of the 4 nearest sites on the lattice (eZ) 2 , 
with probability 1/4, independent of all other particles. Then all particles which 
moved outside D £ are killed. If a site in D £ is occupied by particles of several 
types, then two particles of different types are chosen randomly and also are killed. 
We repeat the procedure, killing pairs of particles of different type occupying the 
same site until there are no sites in D e with more than one type of particle. Killed 
particles will be replaced with new offspring as follows. For every /c, we choose 
rik (to be defined below) particles of type randomly from among the surviving 
ones and each of these particles splits into two offspring of the same type which 
then occupy the same site as the parent particle. Now we define n^. Let n\ be the 
number of particles of type which died because they moved outside D e . Let n 2 . 
be the number of the pairs of particles which were killed inside D s such that the 
types and lifetimes of the particles involved were 7\) and T 2 ) and Ti > T 2 
(i.e., the particle with type had a shorter lifetime). Let n| be defined just as n\ 
except that we replace the condition T\ > T 2 with the condition 7\ = T 2 . Then we 
set nk = n\ + 2n\ + n\. Note that the total number of particles in our model is 
constant between generations but the number of particles of type Ck can vary, for 
each k. 
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Again, we consider the high density limit distribution for the system. Fix some 
open connected set D C R 2 , let D e = D n (eZ) 2 and assume that at time t = 1 
all particles occupy sites in D £ . Recall that we have the total of j particles which 
belong to m different types Ck- Let the measure ' J ' e of an open subset A of R 2 
be equal to the number of particles of type Ck which are in A at time [s]. 

Fix m > 2 and D C R 2 and let j — > oo, e — > and s — > oo. In the limit, 
for every /c, the measure X^'^'^dx, dy)/j will converge to Ckfk(x,y)dxdy (in other 
words, X^' 2 e (A) /j — > Ckfk{x, y)dxdy for every open set A C R 2 ) where < < 
oo and /& is the first eigenfunction of the Laplacian with zero boundary values on a 
subdomain Dk of D. Because of the dynamics, particles of different types become 
segregated so the subdomains Dk are disjoint and their union is D. 

Our transformation rules have been chosen so that the average lifetimes of 
particles of different types are equal in the limit. For if at a certain time the average 
lifetime of particles of type Ck is smaller than that for type C n , the collisions of the 
particles of these two types will result in an increase of the number of particles of 
type Ck- This will imply the growth of the subregion Dk occupied by particles of 
type Ck and hence their average lifetime will increase. The opposite will be true 
for the particles of type C n and so in the limit the average lifetimes of all types of 
particles will be the same. 

The average lifetime of a particle of type Ck is equal to the inverse of the first 
eigenvalue in Dk- Hence, the first eigenvalue for the Laplacian with zero boundary 
conditions in Dk is the same for every k, in the limit. 

Let (x, y) be a point on the boundary between between two subregions Dk and 
D n and let N be the normal unit vector to the boundary at (x, y) pointing inside 
D k - Note that the normal unit vector N at (x, y) pointing inside D n is the same as 
— N. Then we must have dckfk/dN = —dc n f n /d(N) because the particles of both 
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types are killed on the boundary at the same rate. 

V Limit distribution and Laplacian eigenfunctions 

Let F(x,y)dxdy = F m (x,y)dxdy be equal to the limit of X^j^ (dx, dy)/ j on 
Dfc. In other words, F(x,y) = Ckfk(x,y) on Dk and the constants Ck are such 
that dckfk/dN = —dc n f n /d(N) on the boundary between Dk and _D n , where 
N is the inward normal vector on the boundary of Dk and N = —N. Hence, 
dF/dN = —dF/d(N) on the boundary between Dk and D n . 

Suppose that g is an eigenfunction for the Laplacian in D with zero boundary 
values. The lines where g is equal to zero are called the "nodal lines" and they 
divide D into a number of subregions D k . The function g is differentiable, so we 
must have d\g\/dN = —d\g\/d(N) on the boundary between Dk and D n . Moreover, 
\g\ is the first eigenfunction for the Laplacian on every subregion Dk- This suggests 
that F m may be equal to \g\ for some eigenfunction g of the Laplacian in D. 

A simple example shows that for some D and m, the limit distribution F m 
cannot be equal to \g\ for any eigenfunction g in D. This is the case when an 
odd number of "nodal lines" for F m meet at a single point. The number of nodal 
lines meeting at one point must be even for an eigenfunction since the sign of the 
eigenfunction in adjacent regions defined by its nodal lines must alternate. There 
would be no consistent way of assigning signs to adjacent regions if an odd number 
of them met at an intersection point of nodal lines. Fig. lb illustrates a limit 
distribution for a system with 3 particle types. In this case, there are three nodal 
lines for F m which meet at one point and consequently F m cannot be equal to \g\ 
in this case. 

One may ask, then, when the limit distribution F m for a multi-type particle 
system corresponds to a higher eigenfunction. We concentrated our efforts on one 
particular class of domains, namely rectangles D because in this case, the eigenval- 
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ues and the corresponding eigenfunctions can be explicitly calculated. 

Let D = {(x,y) G R 2 : < x < a, < y < b}. Then all eigenvalues of 
the Laplacian in D with zero boundary values are given by Xj t k = n 2 [(j/ a ) 2 + 
(k/b) 2 }, where j and k are arbitrary integers greater than [15]. The eigenfunction 
corresponding to Xj^ has the form fj^(x,y) = sm((jir/a)x) sin((/c7r/6)y). It may 
happen that A jl) / Cl = Xj 2t k 2 even though j\ ^ j 2 and ki ^ k 2 but this is possible 
for only a countable number of side ratios r = b/a. We can also write Xj^ as 
(7v/a) 2 (j 2 + (k/r) 2 ). 

It is intuitively clear that when the number of types of particles m is constant 
but the side ratio of D is very large then particles of different types will occupy 
m rectangles arranged in a linear order (see for example Fig. la). We will call this 
arrangement "elementary." It corresponds to an eigenfunction fj^ of the Lapla- 
cian with either j = 1 or k = 1. This effect is due to the tendency of different 
populations to separate and an elementary configuration seems to be a natural way 
to achieve maximum separation. It is not so clear what happens when the side 
ratio is moderate. When m is fixed, say m = 3, and the side ratio is close to 1, 
we obtain in computer simulations a configuration illustrated in Fig. lb which does 
not correspond to any eigenfunction. We determined by simulation the critical side 
ratio at which we observe the transition between the elementary configuration and 
a configuration which does not correspond to any eigenfunction. 

VI Computer simulations 

Further discussion of the limiting distributions and eigenvalues will be illus- 
trated by computer simulations so we make a digression to explain our figures. In 
all simulations we took D to be a rectangle. The figures show the regions D and 
the boundaries between the subregions occupied by different particle types. All 
simulations were done for rectangles D with sides b = 100 and 100 < a < 300. 
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Because of memory constraints, the results of the simulations were compressed in 
the following way. Every region D was divided into a number of small identical 
rectangles, usually with side lengths between 5 and 10. The numbers of particles of 
different type were found in every small rectangle and the rectangle was declared of 
type Ck if the number of particles of this type was the greatest of all particle types. 
Only rectangles close to the boundaries between D k } s contained different particles 
types. In our simulations, almost all other rectangles contained only one type of 
particles. 

We have simulated long time behavior of the system in rectangles of different 
side ratios with 100,000 particles. Most simulations ran for 150,000 or 200,000 
timesteps. The starting configurations included "elementary configurations," other 
configurations with polygonal separating lines and totally random configurations. 
We used various initial proportions of different particle types. We did simulations 
with m = 3, 4 and 5 particle types. In each case we determined the critical side ratio 
r m = a/6 at which we observed a transformation of the stationary configuration 
from the elementary configuration to a configuration which did not correspond to 
an eigenfunction. The simulations were peformed in 20 different rectangles. Due 
to time consuming nature of the simulations, the number of independent samples 
varied from 1 to 5 per rectangle. The final configurations for the segregation phases 
were unique and did not depend on the initial configuration except when the side 
ratios were close to the critical values discussed below. 

When the number of particle types is m = 3, the critical side ratio is 1.64 ±0.01 
(Fig. 1). The simulations starting from various initial distributions show that the 
limit distribution is elementary for the ratio 1.65 and it is not for the ratio 1.62. 
In the case of side length ratios 1.63 and 1.64, the particle configuration had a 
tendency to preserve its initial shape if the initial shape was as in Fig. la-b. 
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The results of the simulations are most clear in the case of 4 particle types. 
Each of the simulations was started from an asymmetric configuration. The critical 
ratio is 2.26 ± 0.01. The particle distributions are given in Fig. 2. 

Simulations with 5 particle types (Fig. 3) were also started from asymmetric 
distributions. In this case, the critical side ratio is 2.85 ± 0.01. 

An "asymmetric" initial configuration is illustrated in Fig. 4. 

VII Configurational transition and Laplacian eigenvalues 

We will argue that the critical side ratios obtained from the computer simu- 
lations match exceptionally well the critical rectangle side ratios for the following 
problem. When is it true that the elementary configuration with m subregions 
corresponds to the m-th eigenfunction? We order the eigenfunctions according 
to the their eigenvalues, i.e., the m-th eigenfunction corresponds to m-th smallest 
eigenvalue. 

Recall the formulae for the eigenvalues of the Laplacian given in Sect. 3. We 
have Xj t k = (K/a) 2 (j 2 + (k/r) 2 ) for a rectangle with sides equal to a and b and 
side ratio r = b/a. The elementary configuration is defined by the eigenfuction 
corresponding to Ai >m . Whether Ai, m is the m-th eigenvalue depends only on r (it 
does not otherwise depend on the values of a and b). Note that Ai^ < Ai )m for 
fc<mso Xi,m is the m-th eigenvalue if and only if 

Ai, m < A 2 ,i. (7.1) 

This is equivalent to (section V) 

l 2 + (m/r) 2 < 2 2 + (1/r) 2 . (7.2) 

We take m = 3,4,5 and solve this equation for r to obtain the following critical 
side ratios r m : r 3 = ^%J?> « 1.63, r 4 = w 2.24, r 5 = 2 3 / 2 w 2.83. 
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Since our simulations were done on a discrete lattice, the critical side ratio 
values calculated for the rectangle D in R 2 are only approximate. Eigenfunctions 
for the discrete Laplacian on a rectangle D = {(x, y) G Z 2 : 1 < x < a, 1 < y < b} 
are given by f(x 7 y) = g(x)h(y) where g and h satisfy g(0) = g(a + 1) = 0, h(0) = 
h(b + 1) = 0, and 

g(x - 1) - 2g(x) + g(x + 1) = -X x g(x), 1 < x < a, 
h(y - 1) - 2h(y) + h(y + 1) = -X y h(y), 1 < y < b. 

Then A = X x + X y is the eigenvalue corresponding to the eigenfunction f(x,y) = 
g(x)h(y). If g changes sign j — 1 times and h changes the sign k — 1 times then 
A = Xj : k is a discrete analog of A^. We have the following explicit formulae for the 
eigenfunctions and the corresponding eigenvalues. 

gj{x) = sm(jnx/(a + 1)), 

hk(y) = sin(kiry/(b+l)), 

AJ = 2(l-cos(j7r/(a+l))), 

X y k = 2(1 - cos(kn/ (b +1))). 
The discrete analogues of inequalities (7.1) and (7.2) are 

A m ,i < Ai >2 

and 

cos(m7r/(a + 1)) + cos(tt/(6 + 1)) > cos(yr/(a + 1)) + cos(2tt/(6 + 1)). 

In the case b = 100, the critical values for a in the last inequality are in the following 
intervals, 

163 < a < 164, m = 3, 
224 < a < 225, m = 4, 
284 < a < 285, m = 5. 
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These values match very well the critical side lengths discussed in the previous 
section. 

It is quite intriguing that the configuratinal transition takes place for side ratio 
related to the eigenvalue of the Laplacian (Eq(7.1-2)). It would be interesting to 
find an explanation for this phenomenon. We hope that our results will be useful 
in the future studies of the population dynamics. 
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Figure captions 

Figure 1. Nodal lines for stationary distribution of particles with 3 particle 
types. Each region, separated by solid lines is occupied by only one type of parti- 
cles, (a) The side ratio r% = a/b = 1.64; elementary configuration corresponding 
to the third Laplacian eigenfunction (b) rs = 1.63; configuration close to the tran- 
sition point (non-elementary configuration) (c) = 1; configuration far from the 
transition point (non-elementary configuration). 

Figure 2. Nodal lines for stationary distribution of particles with 4 particle 
types. Each region, separated by solid lines is occupied by only one type of parti- 
cles, (a) The side ratio = a/b = 2.27; elementary configuration corresponding to 
the fourth Laplacian eigenfunction (b) = 2.24; configuration close to the tran- 
sition point (non-elementary configuration) (c) r± = 1; configuration far from the 
transition point (non-elementary configuration). 

Figure 3. Nodal lines for stationary distribution of particles with 5 particle 
types. Each region, separated by solid lines is occupied by only one type of parti- 
cles, (a) The side ratio r$ = a/b = 2.88; elementary configuration corresponding 
to the fifth Laplacian eigenfunction (b) rs = 2.84; configuration close to the tran- 
sition point (non-elementary configuration) (c) r$ = 1; configuration far from the 
transition point (non-elementary configuration). 
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Figure 4. An "asymmetric" initial configuration with 4 particle types. Con- 
figurations of this type were used as initial configurations for many simulations. 
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